! ******************************************************
! MPI Communication
! ******************************************************

module Communication
    use Input
    use Boundary
    implicit none

contains

    subroutine exchangeIons()
        implicit none
        type(PseudoParticle), pointer :: p, q

        if(ppx==1)then
        else if(ppi==0)then
            call exchangeIonsRight()
 		else if(ppi==ppx-1)then
            call exchangeIonsLeft()
 		else if(MOD(ppi,2)==1)then
            call exchangeIonsLeft()
            call exchangeIonsRight()
		else
            call exchangeIonsRight()
            call exchangeIonsLeft()
 		end if

        if(ppy==1)then
        else if(ppj==0)then
            call exchangeIonsUp()
 		else if(ppj==ppy-1)then
            call exchangeIonsDown()
 		else if(MOD(ppj,2)==1)then
            call exchangeIonsDown()
            call exchangeIonsUp()
		else
            call exchangeIonsUp()
            call exchangeIonsDown()
 		end if
    end subroutine exchangeIons

    subroutine exchangeIonsRight()
        implicit none
        integer(kind=4) :: obj

        obj = ppx*ppj+ppi+1
        call sendIons(obj,MPIListHeadRight)
        call recvIons(obj,IonListHead)
    end subroutine exchangeIonsRight

    subroutine exchangeIonsLeft()
        implicit none
        integer(kind=4) :: obj

        obj = ppx*ppj+ppi-1
        call recvIons(obj,IonListHead)
        call sendIons(obj,MPIListHeadLeft)
    end subroutine exchangeIonsLeft

    subroutine exchangeIonsUp()
        implicit none
        integer(kind=4) :: obj

        obj = ppx*(ppj+1)+ppi
        call sendIons(obj,MPIListHeadUp)
        call recvIons(obj,IonListHead)
    end subroutine exchangeIonsUp

    subroutine exchangeIonsDown()
        implicit none
        integer(kind=4) :: obj

        obj = ppx*(ppj-1)+ppi
        call recvIons(obj,IonListHead)
        call sendIons(obj,MPIListHeadDown)
    end subroutine exchangeIonsDown

    subroutine exchangeDistribution()
        implicit none

        if(ppx==1)then
        else if(ppi==0)then
            call exchangeGridsRight()
 		else if(ppi==ppx-1)then
            call exchangeGridsLeft()
 		else if(MOD(ppi,2)==1)then
            call exchangeGridsLeft()
            call exchangeGridsRight()
		else
            call exchangeGridsRight()
            call exchangeGridsLeft()
 		end if

        if(ppy==1)then
        else if(ppj==0)then
            call exchangeGridsUp()
 		else if(ppj==ppy-1)then
            call exchangeGridsDown()
 		else if(MOD(ppj,2)==1)then
            call exchangeGridsDown()
            call exchangeGridsUp()
		else
            call exchangeGridsUp()
            call exchangeGridsDown()
 		end if
    end subroutine exchangeDistribution

    subroutine exchangeGridsRight()
        implicit none
        call exchangeGridsRightfor(ni)
        call exchangeGridsRightfor(ne)
        call exchangeGridsRightfor(ux)
        call exchangeGridsRightfor(uy)
        call exchangeGridsRightfor(uz)

    end subroutine exchangeGridsRight

    subroutine exchangeGridsLeft()
        implicit none
        call exchangeGridsLeftfor(ni)
        call exchangeGridsLeftfor(ne)
        call exchangeGridsLeftfor(ux)
        call exchangeGridsLeftfor(uy)
        call exchangeGridsLeftfor(uz)

    end subroutine exchangeGridsLeft

    subroutine exchangeGridsUp()
        implicit none
        call exchangeGridsUpfor(ni)
        call exchangeGridsUpfor(ne)
        call exchangeGridsUpfor(ux)
        call exchangeGridsUpfor(uy)
        call exchangeGridsUpfor(uz)

    end subroutine exchangeGridsUp

    subroutine exchangeGridsDown()
        implicit none
        call exchangeGridsDownfor(ni)
        call exchangeGridsDownfor(ne)
        call exchangeGridsDownfor(ux)
        call exchangeGridsDownfor(uy)
        call exchangeGridsDownfor(uz)

    end subroutine exchangeGridsDown

    subroutine exchangeGridsRightfor(G)
        implicit none
        real(kind=8) :: G(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        integer(kind=4) :: obj
        integer(kind=4) :: num, size
        real(kind=8) :: Temp(2+ncx,-1-ncy:Ny+ncy+1)

        obj = ppx*ppj+ppi+1
        num = Ny+3+2*ncy
        size = 2+ncx
        Temp(1:2+ncx,:) = G(Nx-1-ncx:Nx,:)
        call MPI_SEND(Temp,size*num,mpi_real8,obj,1,mpi_comm_world,ierr)
        call MPI_RECV(Temp,size*num,mpi_real8,obj,1,mpi_comm_world,status,ierr)
        G(Nx:Nx+ncx+1,:) = Temp(1:2+ncx,:)
    end subroutine exchangeGridsRightfor

    subroutine exchangeGridsLeftfor(G)
        implicit none
        real(kind=8) :: G(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        integer(kind=4) :: obj
        integer(kind=4) :: num, size
        real(kind=8) :: Temp(2+ncx,-1-ncy:Ny+ncy+1)

        obj = ppx*ppj+ppi-1
        num = Ny+3+2*ncy
        size = 2+ncx
        call MPI_RECV(Temp,size*num,mpi_real8,obj,1,mpi_comm_world,status,ierr)
        G(0,:) = G(0,:)+Temp(2+ncx,:)
        G(-1-ncx:-1,:) = Temp(1:1+ncx,:)

        Temp(1:2+ncx,:) = G(0:ncx+1,:)
        call MPI_SEND(Temp,size*num,mpi_real8,obj,1,mpi_comm_world,ierr)
    end subroutine exchangeGridsLeftfor

    subroutine exchangeGridsUpfor(G)
        implicit none
        real(kind=8) :: G(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        integer(kind=4) :: obj
        integer(kind=4) :: num, size
        real(kind=8) :: Temp(-1-ncx:Nx+ncx+1,2+ncy)

        obj = ppx*(ppj+1)+ppi
        num = Nx+3+2*ncx
        size = 2+ncy
        Temp(:,1:2+ncy) = G(:,Ny-1-ncy:Ny)
        
        call MPI_SEND(Temp,num*size,mpi_real8,obj,1,mpi_comm_world,ierr)
        call MPI_RECV(Temp,num*size,mpi_real8,obj,1,mpi_comm_world,status,ierr)
        G(:,Ny:Ny+ncy+1) = Temp(:,1:2+ncy)
    end subroutine exchangeGridsUpfor

    subroutine exchangeGridsDownfor(G)
        implicit none
        real(kind=8) :: G(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        integer(kind=4) :: obj
        integer(kind=4) :: num, size
        real(kind=8) :: Temp(-1-ncx:Nx+ncx+1,2+ncy)

        obj = ppx*(ppj-1)+ppi
        num = Nx+3+2*ncx
        size = 2+ncy
        call MPI_RECV(Temp,num*size,mpi_real8,obj,1,mpi_comm_world,status,ierr)
        G(:,0) = G(:,0)+Temp(:,2+ncy)
        G(:,-1-ncy:-1) = Temp(:,1:1+ncy)

        Temp(:,1:2+ncy) = G(:,0:ncy+1)
        call MPI_SEND(Temp,num*size,mpi_real8,obj,1,mpi_comm_world,ierr)
    end subroutine exchangeGridsDownfor

    subroutine smoothGridOnly(U)
        implicit none
        real(kind=8) :: U(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: G(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        integer(kind=4) :: i, j, istart, jstart
        real(kind=8) :: Temp(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        istart = 0
        jstart = 0
        G = U
        if(loc(U)==loc(Bx) .or. loc(U)==loc(Bfx) .or. loc(U)==loc(Bwx) &
        .or. loc(U)==loc(Ey)) then
            istart = 0
            jstart = 1
            G(-1-ncx:Nx+ncx+1,-ncy:Ny+ncy+1) = getBxGrid(U)
        elseif(loc(U)==loc(By) .or. loc(U)==loc(Bfy) .or. loc(U)==loc(Bwy) &
            .or. loc(U)==loc(Ex)) then
            istart = 1
            jstart = 0
            G(-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1) = getByGrid(U)
        elseif(loc(U)==loc(Bz) .or. loc(U)==loc(Bfz) .or. loc(U)==loc(Bwz)) then
            istart = 1
            jstart = 1
            G(-ncx:Nx+ncx+1,-ncy:Ny+ncy+1) = getBzGrid(U)
        
        end if
        
        Temp = 0.0d0
        do j=jstart,Ny
            do i=istart,Nx
                Temp(i,j) = G(i,j)/4.0d0+&
                        (G(i-1,j)+G(i+1,j)+G(i,j-1)+G(i,j+1))/8.0d0+&
                        (G(i-1,j-1)+G(i+1,j-1)+G(i-1,j+1)+G(i+1,j+1))/16.0d0
                ! Temp(i,j) = (G(i,j)+&
                ! G(i-1,j)+G(i+1,j)+G(i,j-1)+G(i,j+1)+&
                ! G(i-1,j-1)+G(i+1,j-1)+G(i-1,j+1)+G(i+1,j+1))/9.0d0
            end do
        end do
        
        if(loc(U)==loc(Bx) .or. loc(U)==loc(Bfx) .or. loc(U)==loc(Bwx) &
        .or. loc(U)==loc(Ey)) then
            U(0:Nx,1:Ny) = Temp(0:Nx,1:Ny)
        elseif(loc(U)==loc(By) .or. loc(U)==loc(Bfy) .or. loc(U)==loc(Bwy) &
            .or. loc(U)==loc(Ex)) then
            U(1:Nx,0:Ny) = Temp(1:Nx,0:Ny)
        elseif(loc(U)==loc(Bz) .or. loc(U)==loc(Bfz) .or. loc(U)==loc(Bwz)) then
            U(1:Nx,1:Ny) = Temp(1:Nx,1:Ny)
        else
            U(0:Nx,0:Ny) = Temp(0:Nx,0:Ny)
        end if
    
    end subroutine smoothGridOnly

    subroutine smoothGrid(U)
        implicit none
        real(kind=8) :: U(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        integer(kind=4) :: i, j
        real(kind=8) :: Temp(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        ! call exchangeGrid(U)
        if(UseConv)then
        else
            call smoothGridOnly(U)
        end if
        call exchangeGrid(U)
    end subroutine smoothGrid

    subroutine exchangeGrid(U)
        implicit none
        real(kind=8) :: U(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)

        call applyGridBoundary_in_exchange(U)
        call applyGridCommunication(U)  ! Communicate changes on the Boudary
    end subroutine exchangeGrid

    subroutine applyGridCommunication(U)
        implicit none
        real(kind=8) :: U(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        integer(kind=4) :: obj
        integer(kind=4) :: num, size
        real(kind=8) :: Tempx(1+ncx,-1-ncy:Ny+ncy+1), Tempy(-1-ncx:Nx+ncx+1,1+ncy)
        size = 1+ncx
        if(ppx==1)then
        else if(ppi==0)then
            obj = ppx*ppj+ppi+1
            num = Ny+3+2*ncy
            if(loc(U)==loc(By) .or. loc(U)==loc(Bfy) .or. loc(U)==loc(Bwy) &
            .or. loc(U)==loc(Bz) .or. loc(U)==loc(Bfz) .or. loc(U)==loc(Bwz) &
            .or. loc(U)==loc(Ex))then
                Tempx = U(Nx-ncx:Nx,:)  ! some Bfield values are on the half Grid
            else
                Tempx = U(Nx-1-ncx:Nx-1,:)
            end if
            call MPI_SEND(Tempx,size*num,mpi_real8,obj,1,mpi_comm_world,ierr)
            call MPI_RECV(Tempx,size*num,mpi_real8,obj,1,mpi_comm_world,status,ierr)
            U(Nx+1:Nx+1+ncx,:) = Tempx
 		else if(ppi==ppx-1)then
            obj = ppx*ppj+ppi-1
            num = Ny+3+2*ncy
            call MPI_RECV(Tempx,size*num,mpi_real8,obj,1,mpi_comm_world,status,ierr)
            U(-1-ncx:-1,:) = Tempx
            Tempx = U(1:1+ncx,:)
            call MPI_SEND(Tempx,size*num,mpi_real8,obj,1,mpi_comm_world,ierr)
 		else if(MOD(ppi,2)==1)then
            obj = ppx*ppj+ppi-1
            num = Ny+3+2*ncy
            call MPI_RECV(Tempx,size*num,mpi_real8,obj,1,mpi_comm_world,status,ierr)
            U(-1-ncx:-1,:) = Tempx
            Tempx = U(1:1+ncx,:)
            call MPI_SEND(Tempx,size*num,mpi_real8,obj,1,mpi_comm_world,ierr)

            obj = ppx*ppj+ppi+1
            if(loc(U)==loc(By) .or. loc(U)==loc(Bfy) .or. loc(U)==loc(Bwy) &
            .or. loc(U)==loc(Bz) .or. loc(U)==loc(Bfz) .or. loc(U)==loc(Bwz) &
            .or. loc(U)==loc(Ex))then
                Tempx = U(Nx-ncx:Nx,:)
            else
                Tempx = U(Nx-1-ncx:Nx-1,:)
            end if
            call MPI_SEND(Tempx,size*num,mpi_real8,obj,1,mpi_comm_world,ierr)
            call MPI_RECV(Tempx,size*num,mpi_real8,obj,1,mpi_comm_world,status,ierr)
            U(Nx+1:Nx+1+ncx,:) = Tempx
		else
            obj = ppx*ppj+ppi+1
            num = Ny+3+2*ncy
            if(loc(U)==loc(By) .or. loc(U)==loc(Bfy) .or. loc(U)==loc(Bwy) &
            .or. loc(U)==loc(Bz) .or. loc(U)==loc(Bfz) .or. loc(U)==loc(Bwz) &
            .or. loc(U)==loc(Ex))then
                Tempx = U(Nx-ncx:Nx,:)
            else
                Tempx = U(Nx-1-ncx:Nx-1,:)
            end if
            call MPI_SEND(Tempx,size*num,mpi_real8,obj,1,mpi_comm_world,ierr)
            call MPI_RECV(Tempx,size*num,mpi_real8,obj,1,mpi_comm_world,status,ierr)
            U(Nx+1:Nx+1+ncx,:) = Tempx

            obj = ppx*ppj+ppi-1
            call MPI_RECV(Tempx,size*num,mpi_real8,obj,1,mpi_comm_world,status,ierr)
            U(-1-ncx:-1,:) = Tempx
            Tempx = U(1:1+ncx,:)
            call MPI_SEND(Tempx,size*num,mpi_real8,obj,1,mpi_comm_world,ierr)
 		end if
        size = 1+ncy
        if(ppy==1)then
        else if(ppj==0)then
            obj = ppx*(ppj+1)+ppi
            num = Nx+3+2*ncx
            if(loc(U)==loc(Bx) .or. loc(U)==loc(Bfx) .or. loc(U)==loc(Bwx) &
            .or. loc(U)==loc(Bz) .or. loc(U)==loc(Bfz) .or. loc(U)==loc(Bwz) &
            .or. loc(U)==loc(Ey))then
                Tempy = U(:,Ny-ncy:Ny)
            else
                Tempy = U(:,Ny-ncy-1:Ny-1)
            end if
            call MPI_SEND(Tempy,size*num,mpi_real8,obj,1,mpi_comm_world,ierr)
            call MPI_RECV(Tempy,size*num,mpi_real8,obj,1,mpi_comm_world,status,ierr)
            U(:,Ny+1:Ny+ncy+1) = Tempy
 		else if(ppj==ppy-1)then
            obj = ppx*(ppj-1)+ppi
            num = Nx+3+2*ncx
            call MPI_RECV(Tempy,size*num,mpi_real8,obj,1,mpi_comm_world,status,ierr)
            U(:,-1-ncy:-1) = Tempy
            Tempy = U(:,1:1+ncy)
            call MPI_SEND(Tempy,size*num,mpi_real8,obj,1,mpi_comm_world,ierr)
 		else if(MOD(ppj,2)==1)then
            obj = ppx*(ppj-1)+ppi
            num = Nx+3+2*ncx
            call MPI_RECV(Tempy,size*num,mpi_real8,obj,1,mpi_comm_world,status,ierr)
            U(:,-1-ncy:-1) = Tempy
            Tempy = U(:,1:1+ncy)
            call MPI_SEND(Tempy,size*num,mpi_real8,obj,1,mpi_comm_world,ierr)

            obj = ppx*(ppj+1)+ppi
            if(loc(U)==loc(Bx) .or. loc(U)==loc(Bfx) .or. loc(U)==loc(Bwx) &
            .or. loc(U)==loc(Bz) .or. loc(U)==loc(Bfz) .or. loc(U)==loc(Bwz) &
            .or. loc(U)==loc(Ey))then
                Tempy = U(:,Ny-ncy:Ny)
            else
                Tempy = U(:,Ny-ncy-1:Ny-1)
            end if
            call MPI_SEND(Tempy,size*num,mpi_real8,obj,1,mpi_comm_world,ierr)
            call MPI_RECV(Tempy,size*num,mpi_real8,obj,1,mpi_comm_world,status,ierr)
            U(:,Ny+1:Ny+ncy+1) = Tempy
		else
            obj = ppx*(ppj+1)+ppi
            num = Nx+3+2*ncx
            if(loc(U)==loc(Bx) .or. loc(U)==loc(Bfx) .or. loc(U)==loc(Bwx) &
            .or. loc(U)==loc(Bz) .or. loc(U)==loc(Bfz) .or. loc(U)==loc(Bwz) &
            .or. loc(U)==loc(Ey))then
                Tempy = U(:,Ny-ncy:Ny)
            else
                Tempy = U(:,Ny-ncy-1:Ny-1)
            end if
            call MPI_SEND(Tempy,size*num,mpi_real8,obj,1,mpi_comm_world,ierr)
            call MPI_RECV(Tempy,size*num,mpi_real8,obj,1,mpi_comm_world,status,ierr)
            U(:,Ny+1:Ny+ncy+1) = Tempy

            obj = ppx*(ppj-1)+ppi
            call MPI_RECV(Tempy,size*num,mpi_real8,obj,1,mpi_comm_world,status,ierr)
            U(:,-1-ncy:-1) = Tempy
            Tempy = U(:,1:1+ncy)
            call MPI_SEND(Tempy,size*num,mpi_real8,obj,1,mpi_comm_world,ierr)
 		end if
    end subroutine applyGridCommunication
end module Communication
